Re-visiting Density Methods

nclust <- function(x0, y0, radius, n) {
              return(runifdisc(n, radius, centre=c(x0, y0)))
            }
x <- rpoispp(lambda =50)
Q = quadratcount(x)
plot(x, main="")
plot(Q, add=TRUE)

Kernel Density Estimates (KDE)

\[ \begin{equation} \hat{f}(x) = \frac{1}{nh_xh_y} \sum_{i=1}^n k\bigg(\frac{{x-x_i}}{h_x},\frac{{y-y_i}}{h_y} \bigg) \end{equation} \]

Kernel Density Estimates (KDE)

x <- rpoispp(lambda =500)
K1 <- density(x, bw=2)
K2 <- density(x, bw=100)
K3 <- density(x, bw=2, kernel="disc")

Warning in density.ppp(x, bw = 2, kernel = "disc"): Bandwidth selection will be
based on Gaussian kernel

par(mfrow=c(1,3))
plot(K1, main="Bandwidth=2")
plot(K2, main="Bandwidth=100")
plot(K3, main="Bandwidth=2, kernel=disc")

par(mfrow=c(1,1))

Kernel Density Estimates (KDE)

• What if we imagine intensity is a function of another covariate

load(url("https://github.com/mgimond/Spatial/raw/main/Data/ppa.RData"))
ma <- readRDS("file/ma.rds")


# Load an MA.shp polygon shapefile 
    s  <- as(ma, "sf")
    w  <- as.owin(s)
    w.km <- rescale(w, 1000) 

    # Load a starbucks.shp point feature shapefile
    s  <- st_read("file/Starbucks.shp", quiet=TRUE)  
    starbucks  <- as.ppp(s)

Warning in as.ppp.sf(s): only first attribute column is used for marks

    marks(starbucks) <- NULL
    Window(starbucks) <- w     
    
    # Load a pop_sqmile.img population density raster layer
    img  <- raster::raster("file/pop_sqmile.img")
    pop  <- as.im(img)
Q <- quadratcount(starbucks, nx= 6, ny=3)
plot(starbucks, pch=20, cols="grey70", main=NULL)  # Plot points
plot(Q, add=TRUE)  # Add quadrat grid    

plot(log(pop))

Adding a Covariate to KDE

• rhohat computes nonparametric intensity \(\rho\) as a function of a covariate

\[ \begin{equation} \lambda(u) = \rho(Z(u)) \end{equation} \]

starbucks.km <- rescale(starbucks, 1000, "km")
ma.km <- rescale(w, 1000, "km")
pop.km    <- rescale(pop, 1000, "km")
pop.lg <- log(pop)
pop.lg.km <- rescale(pop.lg, 1000, "km")
rho <- rhohat(starbucks.km, pop.lg.km,  method="ratio")

Adding a Covariate to KDE

plot(rho, las=1, main=NULL, legendargs=list(cex=0.8, xpd=TRUE, inset=c(1.01, 0) ))

Warning in min(D[scaledlegbox]): no non-missing arguments to min; returning Inf

pred <- predict(rho)
cl   <- interp.colours(c("lightyellow", "orange" ,"red"), 100) # Create color scheme
plot(pred, col=cl, las=1, main=NULL, gamma = 0.25)

Adding a Covariate to KDE

• We can also think more generatively

• Model explicitly as a Poisson Point Process using ppm

\[ \begin{equation} \lambda(u) = \exp^{Int + \beta X} \end{equation} \]

# Create the Poisson point process model
PPM1 <- ppm(starbucks.km ~ pop.lg.km)

Warning: Values of the covariate 'pop.lg.km' were NA or undefined at 0.57% (4
out of 699) of the quadrature points. Occurred while executing: ppm.ppp(Q =
starbucks.km, trend = ~pop.lg.km, data = NULL, interaction = NULL)

# Plot the relationship

Adding a Covariate to KDE

PPM1

Nonstationary Poisson process

Log intensity:  ~pop.lg.km

Fitted trend coefficients:
(Intercept)   pop.lg.km 
 -13.710551    1.279928 

              Estimate       S.E.    CI95.lo    CI95.hi Ztest      Zval
(Intercept) -13.710551 0.46745489 -14.626746 -12.794356   *** -29.33021
pop.lg.km     1.279928 0.05626785   1.169645   1.390211   ***  22.74705
Problem:
 Values of the covariate 'pop.lg.km' were NA or undefined at 0.57% (4 out of 
699) of the quadrature points

plot(effectfun(PPM1, "pop.lg.km", se.fit=TRUE), main=NULL, 
     las=1, legendargs=list(cex=0.8, xpd=TRUE, inset=c(1.01, 0) ))

Warning in min(D[scaledlegbox]): no non-missing arguments to min; returning Inf

Overlays

• Methods for identifying optimal site selection or suitability

• Apply a common scale to diverse or disimilar outputs

Getting Started

Boolean Overlays

Boolean Overlays

• Reclassifying

• Which types of land are appropriate

nlcd <-  rast(system.file("raster/nlcd.tif", package = "spDataLarge"))
plot(nlcd)

Boolean Overlays

• Which types of land are appropriate?

nlcd.segments <- segregate(nlcd)
names(nlcd.segments) <- levels(nlcd)[[1]][-1,2]
plot(nlcd.segments)

Boolean Overlays

• Which types of land are appropriate?

srtm <- rast(system.file("raster/srtm.tif", package = "spDataLarge"))
slope <- terrain(srtm, v = "slope")
par(mfrow=c(1,2))
plot(srtm)
plot(slope)

par(mfrow=c(1,1))

Boolean Overlays

• Make sure data is aligned!

suit.slope <- slope < 10
suit.landcov <- nlcd.segments["Shrubland"]
suit.slope.match <- project(suit.slope, suit.landcov)
suit <- suit.slope.match + suit.landcov

Boolean Overlays

par(mfrow=c(1,3))
plot(suit.slope)
plot(suit.landcov)
plot(suit)

par(mfrow=c(1,1))

Challenges with Boolean Overlays

1. Assume relationships are really Boolean

2. No measurement error

3. Categorical measurements are known exactly

4. Boundaries are well-represented

A more general approach

• Define a favorability metric

\[ \begin{equation} F(\mathbf{s}) = \prod_{M=1}^{m}X_m(\mathbf{s}) \end{equation} \]

Weighted Linear Combinations

\[ \begin{equation} F(\mathbf{s}) = \frac{\sum_{M=1}^{m}w_mX_m}{\sum_{M=1}^{m}w_i} \end{equation} \]